In [19]:
from __future__ import division
import numpy as np
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
Давайте уточним правило трёх сигм. Утверждение: 99.7% вероятностной массы случайной величины X∼N(μ,σ2) лежит в интервале μ±c⋅σ. Чему равно точное значение константы c? Округлите ответ до четырёх знаков после десятичной точки.
In [5]:
from scipy import stats
print('Answer: %.4f' % stats.norm.ppf(1-0.003/2))
В пятилетнем рандомизированном исследовании Гарвардской медицинской школы 11037 испытуемых через день принимали аспирин, а ещё 11034 — плацебо. Исследование было слепым, то есть, испытуемые не знали, что именно они принимают.
За 5 лет инфаркт случился у 104 испытуемых, принимавших аспирин, и у 189 принимавших плацебо.
Оцените, насколько вероятность инфаркта снижается при приёме аспирина. Округлите ответ до четырёх знаков после десятичной точки.
In [22]:
aspirin_data = np.array( [1 if i<104 else 0 for i in range(11037)] )
placebo_data = np.array( [1 if i<189 else 0 for i in range(11034)] )
In [27]:
prob_infarction_aspirin = aspirin_data.sum() / aspirin_data.shape[0]
prob_infarction_placebo = placebo_data.sum() / placebo_data.shape[0]
print('Infarction probability (aspirin): %.4f' % prob_infarction_aspirin)
print('Infarction probability (placebo): %.4f' % prob_infarction_placebo)
print('Infarction probability decrease: %.4f' % (prob_infarction_placebo - prob_infarction_aspirin))
Постройте теперь 95% доверительный интервал для снижения вероятности инфаркта при приёме аспирина. Чему равна его верхняя граница? Округлите ответ до четырёх знаков после десятичной точки.
In [33]:
def proportions_confint_diff_ind(sample1, sample2, alpha=0.05):
z = stats.norm.ppf(1 - alpha / 2.)
p1 = float(sum(sample1)) / len(sample1)
p2 = float(sum(sample2)) / len(sample2)
left_boundary = (p1 - p2) - z * np.sqrt(p1 * (1 - p1)/ len(sample1) + p2 * (1 - p2)/ len(sample2))
right_boundary = (p1 - p2) + z * np.sqrt(p1 * (1 - p1)/ len(sample1) + p2 * (1 - p2)/ len(sample2))
return (left_boundary, right_boundary)
In [36]:
print('Confidence interval: [%.4f, %.4f]' % proportions_confint_diff_ind(placebo_data, aspirin_data))
Продолжим анализировать данные эксперимента Гарвардской медицинской школы.
Для бернуллиевских случайных величин X∼Ber(p) часто вычисляют величину p1−p, которая называется шансами (odds). Чтобы оценить шансы по выборке, вместо p нужно подставить p^. Например, шансы инфаркта в контрольной группе, принимавшей плацебо, можно оценить как
Оцените, во сколько раз понижаются шансы инфаркта при регулярном приёме аспирина. Округлите ответ до четырёх знаков после десятичной точки.
In [37]:
def odds(data):
p = data.sum() / data.shape[0]
return p / (1 - p)
In [40]:
odds_aspirin = odds(aspirin_data)
print('Odds aspirin: %.4f' % odds_aspirin)
odds_placebo = odds(placebo_data)
print('Odds aspirin: %.4f' % odds_placebo)
In [41]:
print('Times decrease of infarction: %.4f' % (odds_placebo / odds_aspirin))
Величина, которую вы оценили в предыдущем вопросе, называется отношением шансов. Постройте для отношения шансов 95% доверительный интервал с помощью бутстрепа. Чему равна его нижняя граница? Округлите ответ до 4 знаков после десятичной точки.
Чтобы получить в точности такой же доверительный интервал, как у нас:
составьте векторы исходов в контрольной и тестовой выборках так, чтобы в начале шли все единицы, а потом все нули;
установите random seed=0;
сделайте по 1000 псевдовыборок из каждой группы пациентов с помощью функции get_bootstrap_samples.
In [43]:
def get_bootstrap_samples(data, n_samples):
indices = np.random.randint(0, len(data), (n_samples, len(data)))
samples = data[indices]
return samples
In [44]:
def stat_intervals(stat, alpha):
boundaries = np.percentile(stat, [100 * alpha / 2., 100 * (1 - alpha / 2.)])
return boundaries
In [48]:
np.random.seed(0)
odds_aspirin_data = np.array(map(odds, get_bootstrap_samples(aspirin_data, 1000)))
odds_placebo_data = np.array(map(odds, get_bootstrap_samples(placebo_data, 1000)))
In [53]:
print('95%% confidence interval for times decrease of infarction: %s' %
str(stat_intervals(odds_placebo_data / odds_aspirin_data, 0.05)))